import numpy as np
import matplotlib.pyplot as plt

Nx = 32
Ny = 32
L = 1 # 求解域长度
dx = 1 / Nx
dy = 1 / Ny

# 网格位置
x = np.linspace(0, L, Nx, endpoint=False)
y = np.linspace(0, Ny * dy, Ny, endpoint=False)
X, Y = np.meshgrid(x, y, indexing='ij')

# print(x)

# for i in range(Nx):
#     print(i*dx)
#     print(x[i])

# 泊松方程等式右边项
f = np.zeros((Nx,Ny))
rsq = np.zeros((Nx,Ny))
for i in range(Nx):
    for j in range(Ny):
        rsq[i,j] = (x[i] - L/2)**2 + (y[j] - L/2)**2
        sigsq = 0.01
        f[i,j] = np.exp(-rsq[i,j] / (2*sigsq))*(rsq[i,j] - 2*sigsq)/(sigsq**2)

plt.figure()
plt.imshow(f, extent=[0, L, 0, L], origin='lower', cmap='viridis')
plt.colorbar(label='f')
plt.show()